home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 007 / bcdfunc.arc / BCDFUN.INC
Encoding:
Text File  |  1985-07-28  |  4.3 KB  |  202 lines
  1. { functions for BCD version of turbo pascal    }
  2. { includes sqrt, log, ln, exp, sin, cos, arctan}
  3. { written by Randall A. Gacek (75026,1042)     }
  4.  
  5. { algorithms adapted from the book:            }
  6. { Software Manual for the Elementary Functions }
  7. { by William J. Cody, jr. & William Waite      }
  8. { Prentice-Hall 1980                           }
  9.  
  10. { WARNING: parts of the code are specific to   }
  11. { the Turbo representation of BCD reals        }
  12.  
  13. function sqrt(x:real):real;
  14.   var
  15.     n,i,m:integer;
  16.     f,y:real;
  17.     v:record case boolean of
  18.       true:(y:real);
  19.       false:(z:array[1..10] of byte)
  20.       end;
  21.   begin
  22.   if x = 0.0 then
  23.     sqrt:=0.0
  24.   else if x < 0.0 then
  25.     halt
  26.   else
  27.     begin
  28.     v.y:=x;
  29.     n:=v.z[1]-63;
  30.     v.z[1]:=63;
  31.     f:=v.y;
  32.     y:=0.580661+f/2.0-0.086462/(f+0.175241);
  33.     for i:=1 to 2 do
  34.       y:=0.5*(y+f/y);
  35.     y:=y+0.5*(f/y-y);
  36.     if odd(n) then
  37.       begin
  38.       y:=y*0.316227766016837933;
  39.       n:=n+1;
  40.       end;
  41.     m:=n div 2;
  42.     v.y:=y;
  43.     v.z[1]:=v.z[1]+m;
  44.     sqrt:=v.y;
  45.     end;
  46.   end; { sqrt }
  47.  
  48. function log(x:real):real;
  49.   const
  50.     c0= 0.316227766016837933;
  51.     a0=-0.260447002405557636E+2;
  52.     a1= 0.554085912041205931E+2;
  53.     a2=-0.392737410203156250E+2;
  54.     a3= 0.103338571514793865E+2;
  55.     a4=-0.741010784161919239E+0;
  56.     b0=-0.899552077881033117E+2;
  57.     b1= 0.245347618868489348E+3;
  58.     b2=-0.244303035341829542E+3;
  59.     b3= 0.107109789115668009E+3;
  60.     b4=-0.193732345832854786E+2;
  61.     c=  0.868588963806503655;
  62.  
  63.   var
  64.     n:integer;
  65.     xn,f,s,w,aw,bw,rs2,rs:real;
  66.     v:record case boolean of
  67.       true:(y:real);
  68.       false:(z:array[1..10] of byte)
  69.       end;
  70.  
  71.   begin
  72.   if x <= 0.0 then
  73.     halt;
  74.   v.y:=x;
  75.   n:=v.z[1]-63;
  76.   v.z[1]:=63;
  77.   f:=v.y;
  78.   if f <= c0 then
  79.     begin
  80.     n:=n-1;
  81.     f:=f*10.0;
  82.     end;
  83.   s:=((f-0.5)-0.5)/(f+1.0);
  84.   w:=sqr(s);
  85.   aw:= (((a4*w+a3)*w+a2)*w+a1)*w+a0;
  86.   bw:=((((w+b4)*w+b3)*w+b2)*w+b1)*w+b0;
  87.   rs2:=w*aw/bw;
  88.   rs:=s*(c+rs2);
  89.   xn:=n;
  90.   log:=xn+rs;
  91.   end; { log }
  92.  
  93. function ln(x:real):real;
  94.   const
  95.     c3=2.30258509299404568;
  96.  
  97.   begin
  98.   ln:=log(x)*c3;
  99.   end;
  100.  
  101. function exp(x:real):real;
  102.   const
  103.     bigx=147.365445951618923;
  104.     smallx=-145.062860858624878;
  105.     eps=5.0e-19;
  106.     onelnsqrt10=0.868588963806503655;
  107.     c1=1.151;
  108.     c2=2.92546497022842009e-4;
  109.     p0=0.333267029226801611e+6;
  110.     p1=0.100974148724273918E+5;
  111.     p2=0.420414268137450315E+2;
  112.     q0=0.666534058453603223E+6;
  113.     q1=0.757393346159883444E+5;
  114.     q2=0.841243584514154545E+3;
  115.     sqrt10=3.16227766016837933;
  116.   var
  117.     n:integer;
  118.     xn,g,z,gpz,qz,rg:real;
  119.     v:record case boolean of
  120.       true:(y:real);
  121.       false:(z:array[1..10] of byte)
  122.       end;
  123.  
  124.   begin
  125.   if x > bigx then
  126.     halt;
  127.   if x < smallx then
  128.     halt;
  129.   if abs(x) < eps then
  130.     exp:=1.0
  131.   else
  132.     begin
  133.     n:=round(x*onelnsqrt10);
  134.     xn:=n;
  135.     g:=(x-xn*c1)-xn*c2;
  136.     z:=sqr(g);
  137.     gpz:=((p2*z+p1)*z+p0)*g;
  138.     qz:= ((z+q2)*z+q1)*z+q0;
  139.     rg:=(0.5+gpz/(qz-gpz))*2.0;
  140.     if odd(n) then
  141.       if n >= 0 then
  142.         rg:=sqrt10*rg
  143.       else
  144.         rg:=rg/sqrt10;
  145.     n:=n div 2;
  146.     v.y:=rg;
  147.     v.z[1]:=v.z[1]+n;
  148.     exp:=v.y;
  149.     end;
  150.   end; { exp }
  151.  
  152. function sincos(x,y,sgn:real):real;
  153.   const
  154.     ymax=3141592654.0;
  155.     onepi=0.318309886183790672;
  156.     c1= 3.141;                  { pi to 22 digits }
  157.     c2= 0.000592653589793238463;
  158.     eps=1.0e-9;
  159.     r1=-0.166666666666666651e+0;
  160.     r2= 0.833333333333316503E-2;
  161.     r3=-0.198412698412018405E-3;
  162.     r4= 0.275573192101527561E-5;
  163.     r5=-0.250521067982745845E-7;
  164.     r6= 0.160589364903715891E-9;
  165.     r7=-0.764291780689104677E-12;
  166.     r8= 0.272047909578888462E-14;
  167.   var
  168.     xn,f,t,g,rg:real;
  169.   begin
  170.   if y >= ymax then
  171.     halt;
  172.   xn:=y*onepi;
  173.   xn:=int(xn+0.5);
  174.   if frac(xn / 2.0) <> 0.0 then
  175.     sgn:=-sgn;
  176.   if abs(x) <> y then { cos wanted }
  177.     xn:=xn-0.5;
  178.   f:=(abs(x)-xn*c1)-xn*c2;
  179.   if abs(f) < eps then
  180.     t:=f
  181.   else
  182.     begin
  183.     g:=sqr(f);
  184.     rg:=(((((((r8*g+r7)*g+r6)*g+r5)*g+r4)*g+r3)*g+r2)*g+r1)*g;
  185.     t:=f+f*rg;
  186.     end;
  187.   sincos:=sgn*t;
  188.   end; { sincos }
  189.  
  190. function sin(x:real):real;
  191.   begin
  192.   if x < 0.0 then
  193.     sin:=sincos(x,-x,-1.0)
  194.   else
  195.     sin:=sincos(x,x,1.0);
  196.   end;
  197.  
  198. function cos(x:real):real;
  199.   begin
  200.   cos:=sincos(x,abs(x)+1.57079632679489662,1.0);
  201.   end;
  202.